import numpy as np
import matplotlib.pyplot as plt
import scipy.io as sio
from sklearn.datasets import load_digits
"""
程序里主要包括了以下几个部分：
1.数据的读入以及特征展示
2.损失函数的定义及优化
3.梯度向量的定义
4.预测函数主体
5.结果展示以及精度评估
"""
data = load_digits()
print(data.keys())

raw_X = data.data#坐标值
raw_y = data.target#标签值
#展示数据形状
print(raw_X.shape)
print(raw_y.shape)

#定义sigmoid函数
def sigmoid(z):
    return 1/(1+np.exp(-z))

def CostFunction(theta,X,y,lamda):
    #theta是一维向量，作矩阵乘法时，计算结果也是一维数组
    R=sigmoid(np.dot(X,theta))
    #X为（1797,465）矩阵，theta为（65,）一维向量，R是(1797,)一维向量
    A=y * np.log(R)
    B=(1-y) * np.log(1-R)
    reg = np.sum(np.power(theta[1:],2)) * (lamda/ (2 * len(X)))
    return -np.sum(A+B)/len(X)+reg #返回损失函数

def gradient_reg(theta,X,y,lamda): #带正则项的梯度向量,防止过度拟合
    A = sigmoid(np.dot(X,theta))
    R = theta[1:]*(lamda / len(X))
    R = np.insert(R,0,values=0,axis=0)
    grad1=(np.dot(X.T,(A-y)))/len(X)
    return grad1+R

X = np.insert(raw_X,0,values=1,axis=1)
print(X.shape)
#损失函数中涉及到y与R的运算，R为（1797，）的一维向量，需要将y也转化为一维向量
y = raw_y.flatten()
print(y.shape)

#用处：Minimization of scalar function of one or more variables.
from scipy.optimize import minimize

def one_vs_all(X,y,lamda,K):  #K是分类器值
    n = X.shape[1]  #获取X的第二个维度的值，此处即列数65
    theta_all = np.zeros((K,n))   #此处为（10,65）矩阵

    for i in range(0,K):#此处K=10，i取1,2,3,...,10
        theta_i = np.zeros(n,)   #以0作为初始值，分别为各分类器运行优化函数
        res = minimize(CostFunction,x0=theta_i,args=(X,y == i,lamda),
                       method='TNC',jac=gradient_reg)  #优化损失函数
        # 将结果保存到theta_all的相应位置
        # 注意i为分类器的标号，存储时需要-1
        theta_all[i-1,:] = res.x
    return theta_all

lamda = 1
K = 10
theta_final = one_vs_all(X,y,lamda,K)
print(theta_final)

def predict(X,theta_final):
    # theta_final是(10,65)，X是(1797,65)，概率最大即为结果
    A = sigmoid(np.dot(X,theta_final.T))
    # np.argmax 用于返回数组中最大数的索引值
    # 一维数组时，返回的是索引值；
    # 二维数组时，返回的是由各索引值对应的一维数组
    # axis=0代表在各行中寻找 axis=1代表在各列中寻找
    k = np.argmax(A,axis=1)
    return k+1

y_predict = predict(X,theta_final)
#因为在分类器中保存的类是1-10，所以这里把标签为10的项更改为0
for i in range(len(y_predict)):
    if y_predict[i]==10:
        y_predict[i]=0
#由于predict方法返回的是分类器的值（索引值加1）组成的一维数组
#所以将返回的预测值和实际标签比较，计算均值就可以得出精度
acc = np.mean(y_predict == y)
print("accuracy: "+str(acc))

#result-plotting function
def showImage(rawX, rawY, k, y_true):
    index = np.random.choice(len(rawX),k*k)
    color = ["Greys_r", "Reds_r"]
    #使用subplots画多图，包含了多个axes
    #sharex sharey代表所有子图共享横轴纵轴
    fig, ax = plt.subplots(k, k, sharex=True, sharey=True, figsize=(k, k))
    #协调每个子图之间的距离
    fig.subplots_adjust(wspace=0.5, hspace=0.5)
    for i in range(k):
        for j in range(k):
            #根据index值找到相应的图片和标签
            image = rawX[index[i * k + j], :]
            label = rawY[index[i * k + j]][0]
            x=0#默认正确的颜色显示是灰色
            if label!=y_true[index[i * k + j]]:
                x=1#错误则显示为红色
            ax[i][j].set_xlabel(label)
            ax[i][j].imshow(image.reshape(8, 8), cmap=color[x])
    plt.xticks([])  # 去掉横坐标值
    plt.yticks([])  # 去掉纵坐标值
    plt.savefig('Classification Result.jpg')
    plt.show()

rawX = data.data
rawy = data.target
y_predict = y_predict.reshape([len(rawy), 1])#y_predict是预测的结果标签
showImage(rawX, y_predict, 8, y)
